!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief sums arrays of real/complex numbers with *much* reduced round-off as compared to
!>      a naive implementation (or the one found in most compiler's SUM intrinsic)
!>      using an implementation of Kahan's algorithm for summing real numbers
!>      that can be used instead of the standard Fortran SUM(array[,mask]).
!>
!>      see also http://en.wikipedia.org/wiki/Kahan_summation_algorithm
!> \note
!>      if the compiler optimises away the 'tricky' bit, no accuracy is gained,
!>      if the compiler uses extended precision inconsistently even worse results might be obtained.
!>      This has not been observed.
!>      This algorithm is not fast, and thus not recommended for cases where round-off is not a
!>      concern but performance is.
!>
!>      the standard intrinsic sum can be 'replaced' using the following use statement
!>
!>      USE kahan_sum, ONLY: sum => kahan_sum
!> \par History
!>      03.2006 [Joost VandeVondele]
!> \author Joost VandeVondele
! **************************************************************************************************
MODULE kahan_sum

   IMPLICIT NONE
   PRIVATE
   PUBLIC :: accurate_dot_product, accurate_sum, accurate_dot_product_2
   INTEGER, PARAMETER :: sp = KIND(0.0), dp = KIND(0.0D0)
   REAL(KIND=sp), PARAMETER :: szero = 0.0_sp
   REAL(KIND=dp), PARAMETER :: dzero = 0.0_dp
   COMPLEX(KIND=sp), PARAMETER :: czero = (0.0_sp, 0.0_sp)
   COMPLEX(KIND=dp), PARAMETER :: zzero = (0.0_dp, 0.0_dp)
   INTEGER, PARAMETER :: dblksize = 8

   INTERFACE accurate_sum
      MODULE PROCEDURE &
         kahan_sum_s1, kahan_sum_d1, kahan_sum_c1, kahan_sum_z1, &
         kahan_sum_s2, kahan_sum_d2, kahan_sum_c2, kahan_sum_z2, &
         kahan_sum_s3, kahan_sum_d3, kahan_sum_c3, kahan_sum_z3, &
         kahan_sum_s4, kahan_sum_d4, kahan_sum_c4, kahan_sum_z4, &
         kahan_sum_s5, kahan_sum_d5, kahan_sum_c5, kahan_sum_z5, &
         kahan_sum_s6, kahan_sum_d6, kahan_sum_c6, kahan_sum_z6, &
         kahan_sum_s7, kahan_sum_d7, kahan_sum_c7, kahan_sum_z7
   END INTERFACE accurate_sum

   INTERFACE accurate_dot_product
      MODULE PROCEDURE &
         kahan_dot_product_d1, &
         kahan_dot_product_s2, kahan_dot_product_d2, kahan_dot_product_z2, &
         kahan_dot_product_d3, kahan_dot_product_masked_d3
   END INTERFACE accurate_dot_product

   INTERFACE accurate_dot_product_2
      MODULE PROCEDURE kahan_blocked_dot_product_d1
   END INTERFACE

CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_s1(array, mask) RESULT(ks)
      REAL(KIND=sp), DIMENSION(:), INTENT(IN)            :: array
      LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL        :: mask
      REAL(KIND=sp)                                      :: ks

      INTEGER                                            :: i1
      REAL(KIND=sp)                                      :: c, t, y

      ks = szero; t = szero; y = szero; c = szero

      IF (PRESENT(mask)) THEN
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1)) THEN
               y = array(i1) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
      ELSE
         DO i1 = 1, SIZE(array, 1)
            y = array(i1) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
      END IF
   END FUNCTION kahan_sum_s1

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_d1(array, mask) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: array
      LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL        :: mask
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i1
      REAL(KIND=dp)                                      :: c, t, y

      ks = dzero; t = dzero; y = dzero; c = dzero

      IF (PRESENT(mask)) THEN
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1)) THEN
               y = array(i1) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
      ELSE
         DO i1 = 1, SIZE(array, 1)
            y = array(i1) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
      END IF
   END FUNCTION kahan_sum_d1

! **************************************************************************************************
!> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
!> \param array1 array of real numbers
!> \param array2 another array of real numbers
!> \return dot product
! **************************************************************************************************
   PURE FUNCTION kahan_dot_product_d1(array1, array2) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:), INTENT(in)            :: array1, array2
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i, n
      REAL(KIND=dp), DIMENSION(dblksize)                 :: c, ks_local, t, y

      t = dzero; y = dzero; c = dzero; ks_local = dzero

      n = SIZE(array1)
      DO i = 1, MOD(n, dblksize)
         y(1) = array1(i)*array2(i) - c(1)
         t(1) = ks_local(1) + y(1)
         c(1) = (t(1) - ks_local(1)) - y(1)
         ks_local(1) = t(1)
      END DO
      DO i = MOD(n, dblksize) + 1, n, dblksize
         y = array1(i:i + (dblksize - 1))*array2(i:i + (dblksize - 1)) - c
         t = ks_local + y
         c = (t - ks_local) - y
         ks_local = t
      END DO
      DO i = 2, dblksize
         y(1) = ks_local(i) - (c(1) + c(i))
         t(1) = ks_local(1) + y(1)
         c(1) = (t(1) - ks_local(1)) - y(1)
         ks_local(1) = t(1)
      END DO
      ks = ks_local(1)
   END FUNCTION kahan_dot_product_d1

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_c1(array, mask) RESULT(ks)
      COMPLEX(KIND=sp), DIMENSION(:), INTENT(IN)         :: array
      LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL        :: mask
      COMPLEX(KIND=sp)                                   :: ks

      COMPLEX(KIND=sp)                                   :: c, t, y
      INTEGER                                            :: i1

      ks = czero; t = czero; y = czero; c = czero

      IF (PRESENT(mask)) THEN
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1)) THEN
               y = array(i1) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
      ELSE
         DO i1 = 1, SIZE(array, 1)
            y = array(i1) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
      END IF
   END FUNCTION kahan_sum_c1

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_z1(array, mask) RESULT(ks)
      COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN)         :: array
      LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL        :: mask
      COMPLEX(KIND=dp)                                   :: ks

      COMPLEX(KIND=dp)                                   :: c, t, y
      INTEGER                                            :: i1

      ks = zzero; t = zzero; y = zzero; c = zzero

      IF (PRESENT(mask)) THEN
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1)) THEN
               y = array(i1) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
      ELSE
         DO i1 = 1, SIZE(array, 1)
            y = array(i1) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
      END IF
   END FUNCTION kahan_sum_z1

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_s2(array, mask) RESULT(ks)
      REAL(KIND=sp), DIMENSION(:, :), INTENT(IN)         :: array
      LOGICAL, DIMENSION(:, :), INTENT(IN), OPTIONAL     :: mask
      REAL(KIND=sp)                                      :: ks

      INTEGER                                            :: i1, i2
      REAL(KIND=sp)                                      :: c, t, y

      ks = szero; t = szero; y = szero; c = szero

      IF (PRESENT(mask)) THEN
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2)) THEN
               y = array(i1, i2) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
      ELSE
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_s2

! **************************************************************************************************
!> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
!> \param array1 2-D array of real numbers
!> \param array2 another 2-D array of real numbers
!> \return dot product
! **************************************************************************************************
   PURE FUNCTION kahan_dot_product_s2(array1, array2) RESULT(ks)
      REAL(KIND=sp), DIMENSION(:, :), INTENT(in)         :: array1, array2
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i1, i2, n1, n2
      REAL(KIND=dp)                                      :: c, t, y

      ks = dzero; t = dzero; y = dzero; c = dzero

      n1 = SIZE(array1, 1)
      n2 = SIZE(array1, 2)
      DO i2 = 1, n2
         DO i1 = 1, n1
            y = REAL(array1(i1, i2), dp)*REAL(array2(i1, i2), dp) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
      END DO
   END FUNCTION kahan_dot_product_s2

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_d2(array, mask) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: array
      LOGICAL, DIMENSION(:, :), INTENT(IN), OPTIONAL     :: mask
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i1, i2
      REAL(KIND=dp)                                      :: c, t, y

      ks = dzero; t = dzero; y = dzero; c = dzero

      IF (PRESENT(mask)) THEN
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2)) THEN
               y = array(i1, i2) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
      ELSE
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_d2

! **************************************************************************************************
!> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
!> \param array1 2-D array of real numbers
!> \param array2 another 2-D array of real numbers
!> \return dot product
! **************************************************************************************************
   PURE FUNCTION kahan_dot_product_d2(array1, array2) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:, :), INTENT(in)         :: array1, array2
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i1, i2, n1, n2
      REAL(KIND=dp)                                      :: c, t, y

      ks = dzero; t = dzero; y = dzero; c = dzero

      n1 = SIZE(array1, 1)
      n2 = SIZE(array1, 2)
      DO i2 = 1, n2
         DO i1 = 1, n1
            y = array1(i1, i2)*array2(i1, i2) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
      END DO
   END FUNCTION kahan_dot_product_d2

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_c2(array, mask) RESULT(ks)
      COMPLEX(KIND=sp), DIMENSION(:, :), INTENT(IN)      :: array
      LOGICAL, DIMENSION(:, :), INTENT(IN), OPTIONAL     :: mask
      COMPLEX(KIND=sp)                                   :: ks

      COMPLEX(KIND=sp)                                   :: c, t, y
      INTEGER                                            :: i1, i2

      ks = czero; t = czero; y = czero; c = czero

      IF (PRESENT(mask)) THEN
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2)) THEN
               y = array(i1, i2) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
      ELSE
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_c2

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_z2(array, mask) RESULT(ks)
      COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN)      :: array
      LOGICAL, DIMENSION(:, :), INTENT(IN), OPTIONAL     :: mask
      COMPLEX(KIND=dp)                                   :: ks

      COMPLEX(KIND=dp)                                   :: c, t, y
      INTEGER                                            :: i1, i2

      ks = zzero; t = zzero; y = zzero; c = zzero

      IF (PRESENT(mask)) THEN
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2)) THEN
               y = array(i1, i2) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
      ELSE
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_z2

! **************************************************************************************************
!> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
!> \param array1 2-D array of complex numbers
!> \param array2 another 2-D array of complex numbers
!> \return dot product
! **************************************************************************************************
   PURE FUNCTION kahan_dot_product_z2(array1, array2) RESULT(ks)
      COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(in)      :: array1, array2
      COMPLEX(KIND=dp)                                   :: ks

      COMPLEX(KIND=dp)                                   :: c, t, y
      INTEGER                                            :: i1, i2, n1, n2

      ks = zzero; t = zzero; y = zzero; c = zzero

      n1 = SIZE(array1, 1)
      n2 = SIZE(array1, 2)
      DO i2 = 1, n2
         DO i1 = 1, n1
            y = array1(i1, i2)*array2(i1, i2) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
      END DO
   END FUNCTION kahan_dot_product_z2

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_s3(array, mask) RESULT(ks)
      REAL(KIND=sp), DIMENSION(:, :, :), INTENT(IN)      :: array
      LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL  :: mask
      REAL(KIND=sp)                                      :: ks

      INTEGER                                            :: i1, i2, i3
      REAL(KIND=sp)                                      :: c, t, y

      ks = szero; t = szero; y = szero; c = szero

      IF (PRESENT(mask)) THEN
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3)) THEN
               y = array(i1, i2, i3) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
      ELSE
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_s3

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_d3(array, mask) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: array
      LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL  :: mask
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i1, i2, i3
      REAL(KIND=dp)                                      :: c, t, y

      ks = dzero; t = dzero; y = dzero; c = dzero

      IF (PRESENT(mask)) THEN
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3)) THEN
               y = array(i1, i2, i3) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
      ELSE
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_d3

! **************************************************************************************************
!> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
!> \param array1 3-D array of real numbers
!> \param array2 another 3-D array of real numbers
!> \return dot product
! **************************************************************************************************
   PURE FUNCTION kahan_dot_product_d3(array1, array2) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(in)      :: array1, array2
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i1, i2, i3, n1, n2, n3
      REAL(KIND=dp)                                      :: c, t, y

      ks = dzero; t = dzero; y = dzero; c = dzero

      n1 = SIZE(array1, 1)
      n2 = SIZE(array1, 2)
      n3 = SIZE(array1, 3)
      DO i3 = 1, n3
         DO i2 = 1, n2
            DO i1 = 1, n1
               y = array1(i1, i2, i3)*array2(i1, i2, i3) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END DO
         END DO
      END DO
   END FUNCTION kahan_dot_product_d3

! **************************************************************************************************
!> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
!>        a mask array determines which product array points to include in the sum
!> \param array1 the first input array to compute the product array
!> \param array2 the second input array to compute the product array
!> \param mask the mask array
!> \param th screening threshold: only array points where the value of mask is greater than th are
!>           included in the sum
!> \return the result of summation
! **************************************************************************************************
   PURE FUNCTION kahan_dot_product_masked_d3(array1, array2, mask, th) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
         POINTER                                         :: array1, array2, mask
      REAL(KIND=dp), INTENT(in)                          :: th
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i1, i2, i3
      REAL(KIND=dp)                                      :: c, t, y

      ks = dzero; t = dzero; y = dzero; c = dzero
      DO i3 = LBOUND(mask, 3), UBOUND(mask, 3)
      DO i2 = LBOUND(mask, 2), UBOUND(mask, 2)
      DO i1 = LBOUND(mask, 1), UBOUND(mask, 1)
         IF (mask(i1, i2, i3) .GT. th) THEN
            y = array1(i1, i2, i3)*array2(i1, i2, i3) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END IF
      END DO
      END DO
      END DO
   END FUNCTION kahan_dot_product_masked_d3

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_c3(array, mask) RESULT(ks)
      COMPLEX(KIND=sp), DIMENSION(:, :, :), INTENT(IN)   :: array
      LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL  :: mask
      COMPLEX(KIND=sp)                                   :: ks

      COMPLEX(KIND=sp)                                   :: c, t, y
      INTEGER                                            :: i1, i2, i3

      ks = czero; t = czero; y = czero; c = czero

      IF (PRESENT(mask)) THEN
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3)) THEN
               y = array(i1, i2, i3) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
      ELSE
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_c3

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_z3(array, mask) RESULT(ks)
      COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN)   :: array
      LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL  :: mask
      COMPLEX(KIND=dp)                                   :: ks

      COMPLEX(KIND=dp)                                   :: c, t, y
      INTEGER                                            :: i1, i2, i3

      ks = zzero; t = zzero; y = zzero; c = zzero

      IF (PRESENT(mask)) THEN
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3)) THEN
               y = array(i1, i2, i3) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
      ELSE
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_z3

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_s4(array, mask) RESULT(ks)
      REAL(KIND=sp), DIMENSION(:, :, :, :), INTENT(IN)   :: array
      LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      REAL(KIND=sp)                                      :: ks

      INTEGER                                            :: i1, i2, i3, i4
      REAL(KIND=sp)                                      :: c, t, y

      ks = szero; t = szero; y = szero; c = szero

      IF (PRESENT(mask)) THEN
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4)) THEN
               y = array(i1, i2, i3, i4) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_s4

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_d4(array, mask) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: array
      LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i1, i2, i3, i4
      REAL(KIND=dp)                                      :: c, t, y

      ks = dzero; t = dzero; y = dzero; c = dzero

      IF (PRESENT(mask)) THEN
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4)) THEN
               y = array(i1, i2, i3, i4) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_d4

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_c4(array, mask) RESULT(ks)
      COMPLEX(KIND=sp), DIMENSION(:, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      COMPLEX(KIND=sp)                                   :: ks

      COMPLEX(KIND=sp)                                   :: c, t, y
      INTEGER                                            :: i1, i2, i3, i4

      ks = czero; t = czero; y = czero; c = czero

      IF (PRESENT(mask)) THEN
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4)) THEN
               y = array(i1, i2, i3, i4) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_c4

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_z4(array, mask) RESULT(ks)
      COMPLEX(KIND=dp), DIMENSION(:, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      COMPLEX(KIND=dp)                                   :: ks

      COMPLEX(KIND=dp)                                   :: c, t, y
      INTEGER                                            :: i1, i2, i3, i4

      ks = zzero; t = zzero; y = zzero; c = zzero

      IF (PRESENT(mask)) THEN
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4)) THEN
               y = array(i1, i2, i3, i4) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_z4

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_s5(array, mask) RESULT(ks)
      REAL(KIND=sp), DIMENSION(:, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      REAL(KIND=sp)                                      :: ks

      INTEGER                                            :: i1, i2, i3, i4, i5
      REAL(KIND=sp)                                      :: c, t, y

      ks = szero; t = szero; y = szero; c = szero

      IF (PRESENT(mask)) THEN
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5)) THEN
               y = array(i1, i2, i3, i4, i5) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_s5

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_d5(array, mask) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i1, i2, i3, i4, i5
      REAL(KIND=dp)                                      :: c, t, y

      ks = dzero; t = dzero; y = dzero; c = dzero

      IF (PRESENT(mask)) THEN
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5)) THEN
               y = array(i1, i2, i3, i4, i5) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_d5

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_c5(array, mask) RESULT(ks)
      COMPLEX(KIND=sp), DIMENSION(:, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      COMPLEX(KIND=sp)                                   :: ks

      COMPLEX(KIND=sp)                                   :: c, t, y
      INTEGER                                            :: i1, i2, i3, i4, i5

      ks = czero; t = czero; y = czero; c = czero

      IF (PRESENT(mask)) THEN
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5)) THEN
               y = array(i1, i2, i3, i4, i5) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_c5

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_z5(array, mask) RESULT(ks)
      COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      COMPLEX(KIND=dp)                                   :: ks

      COMPLEX(KIND=dp)                                   :: c, t, y
      INTEGER                                            :: i1, i2, i3, i4, i5

      ks = zzero; t = zzero; y = zzero; c = zzero

      IF (PRESENT(mask)) THEN
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5)) THEN
               y = array(i1, i2, i3, i4, i5) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_z5

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_s6(array, mask) RESULT(ks)
      REAL(KIND=sp), DIMENSION(:, :, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      REAL(KIND=sp)                                      :: ks

      INTEGER                                            :: i1, i2, i3, i4, i5, i6
      REAL(KIND=sp)                                      :: c, t, y

      ks = szero; t = szero; y = szero; c = szero

      IF (PRESENT(mask)) THEN
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5, i6)) THEN
               y = array(i1, i2, i3, i4, i5, i6) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5, i6) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_s6

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_d6(array, mask) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:, :, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i1, i2, i3, i4, i5, i6
      REAL(KIND=dp)                                      :: c, t, y

      ks = dzero; t = dzero; y = dzero; c = dzero

      IF (PRESENT(mask)) THEN
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5, i6)) THEN
               y = array(i1, i2, i3, i4, i5, i6) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5, i6) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_d6

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_c6(array, mask) RESULT(ks)
      COMPLEX(KIND=sp), DIMENSION(:, :, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      COMPLEX(KIND=sp)                                   :: ks

      COMPLEX(KIND=sp)                                   :: c, t, y
      INTEGER                                            :: i1, i2, i3, i4, i5, i6

      ks = czero; t = czero; y = czero; c = czero

      IF (PRESENT(mask)) THEN
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5, i6)) THEN
               y = array(i1, i2, i3, i4, i5, i6) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5, i6) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_c6

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_z6(array, mask) RESULT(ks)
      COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
         OPTIONAL                                        :: mask
      COMPLEX(KIND=dp)                                   :: ks

      COMPLEX(KIND=dp)                                   :: c, t, y
      INTEGER                                            :: i1, i2, i3, i4, i5, i6

      ks = zzero; t = zzero; y = zzero; c = zzero

      IF (PRESENT(mask)) THEN
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5, i6)) THEN
               y = array(i1, i2, i3, i4, i5, i6) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5, i6) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_z6

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_s7(array, mask) RESULT(ks)
      REAL(KIND=sp), DIMENSION(:, :, :, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
         INTENT(IN), OPTIONAL                            :: mask
      REAL(KIND=sp)                                      :: ks

      INTEGER                                            :: i1, i2, i3, i4, i5, i6, i7
      REAL(KIND=sp)                                      :: c, t, y

      ks = szero; t = szero; y = szero; c = szero

      IF (PRESENT(mask)) THEN
         DO i7 = 1, SIZE(array, 7)
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
               y = array(i1, i2, i3, i4, i5, i6, i7) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i7 = 1, SIZE(array, 7)
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5, i6, i7) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_s7

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_d7(array, mask) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:, :, :, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
         INTENT(IN), OPTIONAL                            :: mask
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: i1, i2, i3, i4, i5, i6, i7
      REAL(KIND=dp)                                      :: c, t, y

      ks = dzero; t = dzero; y = dzero; c = dzero

      IF (PRESENT(mask)) THEN
         DO i7 = 1, SIZE(array, 7)
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
               y = array(i1, i2, i3, i4, i5, i6, i7) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i7 = 1, SIZE(array, 7)
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5, i6, i7) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_d7

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_c7(array, mask) RESULT(ks)
      COMPLEX(KIND=sp), DIMENSION(:, :, :, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
         INTENT(IN), OPTIONAL                            :: mask
      COMPLEX(KIND=sp)                                   :: ks

      COMPLEX(KIND=sp)                                   :: c, t, y
      INTEGER                                            :: i1, i2, i3, i4, i5, i6, i7

      ks = czero; t = czero; y = czero; c = czero

      IF (PRESENT(mask)) THEN
         DO i7 = 1, SIZE(array, 7)
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
               y = array(i1, i2, i3, i4, i5, i6, i7) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i7 = 1, SIZE(array, 7)
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5, i6, i7) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_c7

! **************************************************************************************************
!> \brief ...
!> \param array ...
!> \param mask ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION kahan_sum_z7(array, mask) RESULT(ks)
      COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :, :, :), &
         INTENT(IN)                                      :: array
      LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
         INTENT(IN), OPTIONAL                            :: mask
      COMPLEX(KIND=dp)                                   :: ks

      COMPLEX(KIND=dp)                                   :: c, t, y
      INTEGER                                            :: i1, i2, i3, i4, i5, i6, i7

      ks = zzero; t = zzero; y = zzero; c = zzero

      IF (PRESENT(mask)) THEN
         DO i7 = 1, SIZE(array, 7)
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
               y = array(i1, i2, i3, i4, i5, i6, i7) - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END IF
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      ELSE
         DO i7 = 1, SIZE(array, 7)
         DO i6 = 1, SIZE(array, 6)
         DO i5 = 1, SIZE(array, 5)
         DO i4 = 1, SIZE(array, 4)
         DO i3 = 1, SIZE(array, 3)
         DO i2 = 1, SIZE(array, 2)
         DO i1 = 1, SIZE(array, 1)
            y = array(i1, i2, i3, i4, i5, i6, i7) - c
            t = ks + y
            c = (t - ks) - y
            ks = t
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
         END DO
      END IF
   END FUNCTION kahan_sum_z7

! **************************************************************************************************
!> \brief computes the accurate sum of blocks of regular dot products
!> \param array1 array of real numbers
!> \param array2 another array of real numbers
!> \param blksize ...
!> \return dot product
! **************************************************************************************************
   FUNCTION kahan_blocked_dot_product_d1(array1, array2, blksize) RESULT(ks)
      REAL(KIND=dp), DIMENSION(:), INTENT(in)            :: array1, array2
      INTEGER, INTENT(IN), OPTIONAL                      :: blksize
      REAL(KIND=dp)                                      :: ks

      INTEGER                                            :: my_blksize
      REAL(KIND=dp)                                      :: DDOT

      my_blksize = 32
      IF (PRESENT(blksize)) my_blksize = blksize

      IF (my_blksize <= 1) THEN
         ! The original should be faster
         ks = accurate_dot_product(array1, array2)
      ELSE IF (my_blksize >= SIZE(array1)) THEN
         ! Just use standard dot product from BLAS for performance
         ks = DDOT(SIZE(array1), array1(1), 1, array2(1), 1)
      ELSE
         ks = 0.0_dp
         BLOCK
            INTEGER :: i, n, stripesize
            REAL(KIND=dp)                                      :: c, dotproduct, t, y
            t = dzero; y = dzero; c = dzero

            n = SIZE(array1)
            DO i = 1, n, my_blksize
               ! Remove 1 to save an operation in the length
               stripesize = MIN(my_blksize, n - i + 1)
               ! Perform dot product on the given stripe
               dotproduct = DDOT(stripesize, array1(i), 1, array2(i), 1)
               y = dotproduct - c
               t = ks + y
               c = (t - ks) - y
               ks = t
            END DO
         END BLOCK
      END IF
   END FUNCTION kahan_blocked_dot_product_d1

END MODULE kahan_sum
